clc; clear; close all;

%% 参数范围和时间设置
alpha_list = linspace(2, 3, 10);
q10_list = linspace(0, 1, 10);

dt = 0.01;
T_total = 400;
T_trans = 350;
N_steps = round(T_total / dt);
N_trans = round(T_trans / dt);
delta0 = 1e-7;

N_alpha = length(alpha_list);
N_q10 = length(q10_list);
N_total = N_alpha * N_q10;

ClassMapVec = zeros(N_total, 1);

B = [1, -4, -3.5;
     0, 1, 2;
    -1, -4, 1.5];

fprintf('开始并行计算状态分类...\n');
parfor idx = 1:N_total
    ia = mod(idx-1, N_alpha) + 1;
    iq = floor((idx-1) / N_alpha) + 1;

    alpha = alpha_list(ia);
    q10 = q10_list(iq);

    % 构造 A 矩阵
    A = [-alpha, 0, 0;
          0,     0, 0;
          0,     0, 0];

    % 初始状态
    x0 = [1e-6; 0; 0; q10];

    % 系统函数句柄
    f = @(t,x) mCNN_single(t, x, A, B);

    % 数值积分轨迹
    traj = zeros(length(x0), N_steps + 1);
    traj(:,1) = x0;
    x = x0;
    for i = 1:N_steps
        x = RK4(f, x, dt);
        traj(:, i+1) = x;
    end

    % 稳态轨迹截取
    steady_traj = traj(1, N_trans+1:end);

    % 最大李雅普诺夫指数估计
    LE_max = estimate_max_LE_Jacobian(f, x0, dt, T_total, T_trans, delta0, N_steps, N_trans);

    % 状态分类
    class_num = Classify(steady_traj, 10, LE_max);

    ClassMapVec(idx) = class_num;
end
fprintf('计算完成。\n');

% 重塑为二维矩阵 (q10在行，alpha在列)
ClassMap = reshape(ClassMapVec, [N_q10, N_alpha]);

%% 绘制分类图
figure;
imagesc(alpha_list, q10_list, ClassMap');
set(gca, 'YDir', 'normal');

% 设置colormap，顺序和类别一致
colormap([ ...
    0.2 0.4 1.0;   % P0 定点
    0.2 0.6 1.0;   % P1
    0.4 0.6 0.9;   % P2
    0.6 0.6 1.0;   % P3
    0.6 0.7 1.0;   % P4
    0.7 1.0 0.7;   % P5
    1.0 1.0 0.5;   % P6
    1.0 0.6 0.0;   % MP 多周期
    0.6 0.2 0.2;   % CH 混沌
    0.5 0.5 0.5    % UB 不稳定
    ]);
colorbar('Ticks', 1:10, ...
    'TickLabels', {'P0','P1','P2','P3','P4','P5','P6','MP','CH','UB'});

caxis([0.5 10.5]);

xlabel('\alpha');
ylabel('q_{10}');

grid on;

%% 函数

function LE_max = estimate_max_LE_Jacobian(f, x0, dt, t_total, t_trans, delta0, N_steps, N_trans)
    h = dt;
    steps = N_steps;
    trans_steps = N_trans;
    
    x = x0;          % 主轨迹
    dim = length(x0);
    v = delta0 * [1; zeros(dim-1,1)];  % 固定扰动方向
    
    sum_ln = 0;
    count = 0;
    
    for i = 1:steps
        x = RK4(f, x, h);
        J = Jacobian(f, x, 1e-8);
        v = v + h * (J * v);
        
        if mod(i, 10) == 0
            dist = norm(v);
            if i > trans_steps
                sum_ln = sum_ln + log(dist/delta0);
                count = count + 1;
            end
            v = (v / dist) * delta0;
        end
    end
    
    LE_max = sum_ln / (count * h * 10);
    
    if abs(LE_max) > 10
        LE_max = NaN;
    end
end

function J = Jacobian(f, x, eps)
    n = length(x);
    fx = f(0, x);
    m = length(fx);
    J = zeros(m, n);
    for i = 1:n
        dx = zeros(n,1);
        dx(i) = eps;
        J(:, i) = (f(0, x + dx) - fx) / eps;
    end
end
